/*This script generates results at the PLSS level data found in the appendix;
Table B2
Table D2
Table D3
Table D4
Table D5
Table D6

Figure D1
Figure D3

Notes tables Table D1 and D7 are built in PLSS_Main.do 

*/


cd "$data"
set more off

*Bring in the PLSS data
use PLSS_Turbines_Irrigation_JAERE.dta, replace

*Globals drawn on
global xcontrols elev_mean elev_std soil stream3_dist tmean ppt_mean lat lon 
global xend sharecropland sharedeveloped shareforest powerlines_dist
global xwind ave_wind_class max_wind_class
global SE vce(cl alt_county_id)  

/*Table B2*/

global sumstat share_og share_og_neg20 ave_wind_class max_wind_class CPIS_per cpis_ind2 cpis_neighbors irr_mean irrigated sharecropland turbine turbine_count p_year  elev_mean elev_std soil stream3_dist tmean ppt_mean sharedeveloped shareforest powerlines_dist

quietly reghdfe turbine CPIS_per irrigated $xwind $xcontrols $xend , absorb(alt_state_id) $SE
	gen samp=e(sample)

foreach stat in count mean sd min max {
	estpost tabstat $sumstat if samp==1, columns(statistics) statistics(`stat') 
	estadd matrix r=e(`stat')
	eststo tab_plss_des_`stat'
	}
	
	drop samp
	
esttab tab_plss_des_* using "$results/tables/tabB2_des_plss.csv",  label main(r)  nostar  mtitles("Count" "Mean" "St. Dev" "Min" "Max") nonote noobs replace



/*Table D2*/
	local i = 1	
	foreach fe in alt_state_id alt_county_id alt_plssid fe_town_3 fe_town_4 fe_town_9 {
		
		qui reghdfe turbine CPIS_per irrigated $xwind $xcontrols $xend , absorb(`fe') cluster(cluster_km)
			qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		estadd scalar MDV = `mdv'
		estadd scalar Nfe = e(df_a_initial)
		estadd scalar Nclus =  e(N_clust1)
	
	estimates store turb_SE_cpis_irr_`i'
	
		local i=`i'+1
	
		}
	

esttab turb_SE_cpis_irr* using "$results/tables/tabD2_PLSS_turb_SE_rob_irr.csv",  label  keep(CPIS_per irrigated ave_wind_class max_wind_class sharecropland) order(CPIS_per irrigated sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace



/*Table D3*/
*OLS specifications
	local i = 1	
	foreach fe in alt_state_id alt_county_id alt_plssid {
		
		qui reghdfe turbine CPIS_per irrigated $xwind $xcontrols $xend , absorb(`fe') $SE
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		estadd scalar MDV = `mdv'
		estadd scalar Nfe = e(df_a_initial)
		estadd scalar Nclus =  e(N_clust1)
	
		estimates store turb_logit_irr_`i'

		local i=`i'+1
	
		}

* Logit specifications
foreach fe in alt_state_id alt_county_id alt_plssid {

	qui logit turbine CPIS_per irrigated $xwind $xcontrols $xend i.`fe' if turb_`fe'_ind==1, $SE
		qui sum turbine if e(sample) ==1 
		local mdv = r(mean)
		margins, dydx(CPIS_per) 
		local odd=el(r(b),1,1)
		estadd scalar MDV = `mdv'
		estadd scalar Margins = `odd'
		estadd scalar Nfe = e(k)-17
		estadd scalar Nclus =  e(N_clust)
	
	estimates store turb_logit_irr_`i'
	

	local i=`i'+1
	
	}	

esttab turb_logit_irr_* using "$results/tables/tabD3_PLSS_turb_logit_irr.csv",  label  keep(CPIS_per irrigated ave_wind_class max_wind_class sharecropland) order(CPIS_per irrigated sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 pr2 se(a3) b(a3) scalars("MDV Mean DV" "Margins MEM of CPIS" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

/*Table D4*/
*OLS estimates for reference
	local i = 1	
	foreach fe in alt_state_id alt_county_id alt_plssid {
		
		qui reghdfe turbine_count CPIS_per irrigated $xwind $xcontrols $xend , absorb(`fe') $SE
			qui sum turbine_count if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(df_a_initial)
			estadd scalar Nclus =  e(N_clust1)
	
			estimates store turb_count_`i'

		local i=`i'+1
		}

*OLS, only positive CPIS sections
	foreach fe in alt_state_id alt_county_id alt_plssid {

		qui reghdfe turbine_count CPIS_per irrigated $xwind $xcontrols $xend if cpis_ind2==1, absorb(`fe') $SE
			qui sum turbine_count if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(df_a_initial)
			estadd scalar Nclus =  e(N_clust1)
	
			estimates store turb_count_`i'

		local i=`i'+1
		}

*Tobit models	
local i=1
	foreach fe in alt_state_id alt_county_id alt_plssid {

		
		qui tobit turbine_count CPIS_per irrigated $xwind $xcontrols $xend i.`fe' if turb_`fe'_ind==1 , ll(0) vce(cl alt_county_id)
			qui sum turbine_count if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(k)-18
			estadd scalar Nclus =  e(N_clust)
	
		estimates store tobit_turb_count_`i'
	

		local i=`i'+1
		}
	

esttab turb_count_* tobit_turb_count_* using "$results/tables/tabD4_PLSS_turb_count_irr.csv",   label  keep(CPIS_per irrigated ave_wind_class max_wind_class sharecropland) order(CPIS_per irrigated sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 pr2 se(a3) b(a3) scalars("MDV Mean DV" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

/*Table D5*/
gen share_all=1 // creates a variable to help with the loop for the full sample
 
local i = 1
	foreach fe in alt_state_id alt_county_id alt_plssid {
		foreach samp in share_all share_og  {
			qui reghdfe turbine CPIS_per irrigated $xwind $xcontrols $xend if `samp'>0.5, absorb(`fe') $SE
				qui sum turbine if e(sample) ==1 
				local mdv = r(mean)
				estadd scalar MDV = `mdv'
				estadd scalar Nfe = e(df_a_initial)
				estadd scalar Nclus =  e(N_clust1)
	
			estimates store samp_turb_ind_`i'

			local i=`i'+1
			}
		}
	drop share_all
	
esttab samp_turb_ind_* using "$results/tables/tabD5_PLSS_turb_ind_sample_rob_irr.csv", label  keep(CPIS_per irrigated ave_wind_class max_wind_class sharecropland) order(CPIS_per irrigated sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace

/*Figure D1*/
*Note, the example indicators are added ex post using a adobe acrobat pro, pasting in the elements from the top row of Figure D2, provided in the excel spreadsheet fig3_fig5_figD2_CPISlayouts

histogram CPIS_per if CPIS_bin>0 & CPIS_per>.05, ///
	percent w(.01) xline(.2, lcolor(black) lpattern(dash)) xline(.4, lcolor(black) lpattern(dash)) xline(.6, lcolor(black) lpattern(dash)) xline(.8, lcolor(black) lpattern(dash)) ///
		scheme(sj) graphregion(fcolor(white) style(none) color(gs16)) xscale(range(0.05 1)) xlabel(0.1(.1)1) ylabel( , nogrid) legend(off)  xsize(8.5) ysize(3.5)

						graph export  "$results/figures/figD1_CPIS_share_histogram.pdf", replace

						
/*Figure D3 and Table D6*/
	local i = 1	
	foreach fe in alt_state_id alt_county_id alt_plssid {  
		
		qui reghdfe turbine i.CPIS_bin irrigated $xwind $xcontrols $xend , absorb(`fe') $SE
			qui sum turbine if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(df_a_initial)
			estadd scalar Nclus =  e(N_clust1)
	
		estimates store turb_CPIS_ind_`i'
	
		local i=`i'+1
	
		}
	
	
	****CPIS only
	foreach fe in alt_state_id alt_county_id alt_plssid { 

		qui reghdfe turbine CPIS_per irrigated $xwind $xcontrols $xend if cpis_ind2==1, absorb(`fe') $SE
			qui sum turbine if e(sample) ==1 
			local mdv = r(mean)
			estadd scalar MDV = `mdv'
			estadd scalar Nfe = e(df_a_initial)
			estadd scalar Nclus =  e(N_clust1)
	
			estimates store turb_CPIS_ind_`i'

		local i=`i'+1
		}
	
esttab turb_CPIS_ind_* using "$results/tables/tabD6_PLSS_turb_CPIS_ind_irr.csv",  label  keep(irrigated 1.CPIS_bin 2.CPIS_bin 3.CPIS_bin 4.CPIS_bin 5.CPIS_bin CPIS_per ave_wind_class max_wind_class sharecropland) order(CPIS_per 1.CPIS_bin 2.CPIS_bin 3.CPIS_bin 4.CPIS_bin 5.CPIS_bin irrigated sharecropland ave_wind_class max_wind_class) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("MDV Mean DV" "Nfe N Fixed Effects" "Nclus Standard Error Clusters") replace


**plot of county or plssid fixed effects
preserve
	
		
		clear
		estimates restore  turb_CPIS_ind_2
		local binlist  1 2 3 4 5 
		local numbins = wordcount("`binlist'")
			
		set obs `numbins'
		gen bin =.
			
		gen beta=.
		gen se=.
		
		gen id=_n
		
		forvalue i = 1/`numbins' {
			local bin: word `i' of `binlist'
			replace bin = `bin' if id == `i'
			
			replace beta=_b[`bin'.CPIS_bin] if id == `i'
			
			replace se=_se[`bin'.CPIS_bin] if id == `i'
			}
			
		
		gen plus95=beta+1.96*se
		gen minus95=beta-1.96*se

		gen plus90=beta+1.64*se
		gen minus90=beta-1.64*se
		
		label define cpis_num 0 "No CPIS" 1 "< 1 quarter" 2 "1 quarter" 3 "2 quarters" 4 "3 quarters" 5 "4 quarters"
		label val bin cpis_num

		label var bin "Equivalent complete PLSS quarter-sections with CPIS"
		
	
		
		twoway (scatter beta bin, lc(forest_green) mc(forest_green)) (rcap plus90 minus90 bin, lpattern(solid) lcolor(forest_green) msize(large)) , ///
		scheme(sj) graphregion(fcolor(white) style(none) color(gs16)) ytitle("{&Delta}Pr(Turbine)" ) xlabel(,valuelabel angle(-30)) xscale(range(0.5 5.5)) ylabel( , nogrid) yline(0) legend(off) note("90th% Confidence Interval") xsize(6.5) ysize(3.5)
				graph export  "$results/figures/figD3_plss_CPIS_bin_irr.pdf", replace
		
		restore

***Counter factual calculation
preserve
		estimates restore turb_CPIS_ind_2
			gen samp=1 if _est_turb_CPIS_ind_2==1
			
		predict hat_turb if samp==1
		replace CPIS_bin=0 if CPIS_bin>0
		
		replace irrigated=0 if irrigated==1
		
		predict hat_0_turb if samp==1
		
		
		egen hat_sum=total(hat_turb)
		egen hat_0_sum=total(hat_0_turb)
		
		gen hat_dif=hat_0_sum-hat_sum
		
		gen hat_perc=hat_dif/hat_sum
				
		log using "$results/intext/discussion_plss_turbine_counter", text replace

		sum hat_dif hat_perc hat_sum
		
		log close		
	restore
	
	
/*Addiional numbers referenced in the text*/
/*footnote 28*/
foreach fe in alt_county_id alt_plssid  {
		bysort `fe': egen meancount_`fe'=mean(turbine_count)
		gen `fe'_vary=turbine_count-meancount_`fe'
		bysort `fe': egen `fe'_variation=max(`fe'_vary)
	}
	log using "$results/intext/footnote28_within_variation", text replace

	sum alt_county_id_variation if alt_county_id_variation!=0
	sum alt_plssid_variation if alt_plssid_variation!=0
	log close
	
	drop *_vary *_variation meancount_*


********THIS COMPLETES THE PLSS LEVEL OUTPUT***************
